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^ ■ Abstract 

We follow the contraction and evolution of a typical Jupiter-mass clump created by the disk instability 

5_( ' mechanism, and compute the rate of planetesimal capture during this evolution. We show that such a 

, 

clump has a slow contraction phase lasting ~ 3 x 10^ years. By following the trajectories of planetesimals 
as they pass through the envelope of the protoplanet, we compute the cross-section for planetesimal 
capture at all stages of the protoplanet's evolution. We show that the protoplanet can capture a large 
fraction of the solid material in its feeding zone, which will lead to an enrichment of the protoplanet in 
heavy elements. The exact amount of this enrichment depends upon, but is not very sensitive to the size 
and random speed of the planetesimals. 
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1 Introduction 

There are nearly 200 known extrasolar planetary systems. Many of these systems contain planets of several 
Jupiter masses. For those few systems that transit the stellar disk, the radius can be measured, and most of 
these transiting planets appear to be gas giants (see, e.g. Henry et al. 2000; Udalski et al. 2002; Pont et al. 
2004; Alonso et al. 2004). Currently, there are two main theories for forming such planets. The core accretion 
theory (Pollack et al. 1996) argues that a heavy element core is built up by the accretion of planetesimals. As 
the core grows, its ability to accrete gas from the surrounding disk increases. When the core is sufficiently 
massive, there is a rapid accretion of such gas and a giant planet is formed. The big advantage of this 
model is that the same basic mechanism of planetesimal accretion will form terrestrial planets, Uranus and 
Neptune-sized intermediate planets, and giant planets. The disadvantages are that the core mass required to 
form Jupiter is at the upper end of the mass estimated from interior models (Saumon and Guillot 2004), and 
that the time to reach a Jupiter mass is uncomfortably close to the upper limit estimated for lifetime of the 
gas disk (Haisch, Lada, & Lada 2001). Both these objections can be ameliorated, if not entirely removed, if 
the accreted material does not sink entirely to the core (Pollack et al. 1996) and if the opacities are in fact 
lower than those estimated from interstellar grains (Podolak 2003; Hubickyj et al. 2005). 

A competing mechanism for giant planet formation is a local disk instability (Boss 1998a). This model 
suggests that under the right conditions an instability can form in the protoplanetary disk. This instability 
can lead to the creation of a self- gravitating clump of gas and dust. Such clumps can contract to form giant 
gaseous protoplancts (Boss 1997; Boss 1998b). This model has several advantages, among them, that planets 
can be formed quickly, before the nebular gas dissipates. Calculations by Boss (2000) show that the unstable 
disk can break up into giant gaseous protoplanets in ~ 10^ yrs. The instability occurs over dynamical 
timescales (some tens of orbital periods) , and the formation of clumps (future giant gaseous protoplanets) 
can easily occur within the estimated lifetime of most circumstellar disks (Haisch, Lada, & Lada 2001). 

One of the problems with the disk instability model is that the planets formed by this mechanism 
start with a solar abundance of elements. Observations as well as theoretical models, however, indicate 
that Jupiter's envelope (and Saturn's as well) is enriched with heavy elements (Young 2003; Saumon et al. 
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1995). The theoretical estimate for the mass of heavy elements in Jupiter is ~ 20 — SOM^ (Saumon and 
Guillot 2004). A solar composition planet of Jupiter's mass would be expected to have only ~ 6M^ of 
heavy elements. This means that planets created by the disk instability mechanism must have accreted this 
additional material later, presumably as solids. In this paper we follow the evolution of an isolated clump 
and compute the planetesimal accretion rate as a function of time. The calculation is divided into three 
parts: the computation of the evolution of an isolated clump, the computation of the cross section for the 
two body interaction between the protoplanet and a planetesimal, and the actual calculation of the mass 
that is accreted. The details are presented in the following sections. 



2 Planetary Evolution Code 

In our calculation we assume an isolated non-rotating clump. The initial physical parameters of the body 
were chosen to fit the expected initial conditions after the onset of the gravitational instability (Boss 2002). 
We assume that there are no external influences on the body, such as solar heating or disk shear. We follow 
the evolution of this clump using a stellar evolution code developed by one of us (Kovetz). 

This code, which was originally developed for stellar evolution studies, solves the standard equations of 
stellar evolution: 

dm 3 p' ^ ' 

dm 47rr4' ' 
du d 1 dL 

m+''di-p = '-d^' 

dt ^ dm ^ ^ dm ^ ' 

— - V— (5) 

dm dm' 

where Yj is the number fraction of the j'th species, related to the mass fraction Xj by Yj = Xj/Aj, Fj is the 
particle flux of the j'th species, determined by the corresponding coefficient of diffusion aj, V(r, m, L, p, T, Y) 
is the temperature 'gradient' dlogT/dlogp, determined by the Mixing Length Recipe (MLR), and the 
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remaining symbols are in standard notation. At the center r, the energy fliix L and the Fj's all vanish. The 
star's surface is taken to be the photosphere. Thus, the surface boundary conditions are L = A-kt'^ctsbT^, 
Fj = and up = GmTg/r'^, where Tg is the optical depth of the photosphere, which we take to be unity. 

The foregoing equations are replaced by implicit difference equations, which are then solved numerically. 
Instead of using a fixed grid of mass points, the code determines the mass distribution by requiring the 
function 



to change by a constant increment between any two consecutive mass points. With c a numerical constant 
of order unity, the second-order difference equation <Pf = 0, which has the boundary conditions m = at 
the center and m = M at the surface , ensures equal steps of logp, except near the center, where it imposes 



In the absence of any nuclear (or chemical) change, the rates Rj will all vanish, but the heating source q 
may still be positive, e.g. when (and where) accreted planetesimals are slowed by friction. 

As noted above, convection is treated in accordance with MLR. In order to avoid the difficulties connected 
with sudden, instantaneous mixing in convective zones, the code regards convective mixing as a diffusive 
process in a gas of density p, mean velocity Vc (supplied by MLR) and mean free path 1^ = 1.5Hp. The 
artificial convective diffusion coefficient CTc = {A-kt'^ p^Vclc (the same for all species) of course vanishes outside 
the convective zones. Any real diffusion, or settling, process can then be easily incorporated by adding actual 
diffusion coefficients. 

The equation of state tables were kindly provided by D. Saumon and are based on Saumon et al. (1995), 
supplemented at low pressures by our own equation of state, appropriate for a weakly interacting (Debye 
approximation) mixture of gases. The opacity tables were kindly provided by P. Bodenheimer, based on the 
work of Pollack et al. (1989). They include both gas and grain opacity, the latter being based on the size 
distribution relevant for interstellar grains. 

We were able to find an initial quasi-static model of a Jupiter mass object, similar to the initial clump 
in the model of Boss (2002), based on a preliminary model, also provided by Bodenheimer. Table 1 gives 




(6) 



equal steps of vn?/^. 
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the initial parameters of the starting model. The pressure and temperature profiles for the initial model are 
shown in fig. 1. 



Properties of the initial Model 




Radius (cm) 


7.f4 X fO^^ 


Effective temperature (K) 


26.3 


Central Temperature (K) 


3.5f X 10^ 


Photospheric density (g cm^'^) 


1.64 X 10"" 


Central density (g cm~'^) 


2.53 X 10-^ 


Photospheric pressure(dyne cm~^) 


1.55 X 10^2 


Central pressure(dyne cm~^) 


3.17 X 10^ 



Table 1: Properties of the initial model 

10' r . . , , . . . . r- 




Radius [cm] 



Figure 1: Pressure and temperature as a function of radius for the initial model 

Fig. 2 shows evolution of the protoplanet's radius, central temperature, and central pressure with time 
for the first 3 x lO^yr. At this point the central temperature reaches ~ 2000 K and the molecular hydrogen 
dissociates sufficiently so that the resultant energy sink triggers a rapid collapse of the body. The evolutionary 
track is similar to one calculated by Bodenheimer at al. (1980), for a IMj protoplanet with an initial central 
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density of 10~^g cm"^ and an initial central temperature of 100 K. For comparison our lAIj protoplanet 
initially has an central density of ~ 2 x 10~^g cm~^, a radius of 0.5 AU, and a central temperature of ~ 350 
K. The evolution of the body was computed assuming no planetesimal capture, so that the composition of 
the protoplanet was constant. 
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Figure 2: Evolution with time of the central pressure (solid curve), the central temperature (dashed curve), 
and the radius of the protoplanet in Jupiter radii (dash-dot curve). 



3 Two-Body Capture Cross Section 

Once the structure of the planet is determined, one can, in principle, compute its capture cross-section. 
If this is known, the time-dependent rate of mass increase due to the capture of planetesimals is given by 
(Safronov 1969) 

^ = (7) 

where S](i) is the protoplanet's cross section for planetesimal capture, cr is the surface density of solid 
material, and f2 is the orbital frequency of the protoplanet. 

To compute the capture cross section at each stage of the evolution, we followed the trajectory of a 
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planetesimal as it passed through the protoplanetary atmosphere. At each step of the trajectory we computed 
the motion of the planetesimal in response to the gas drag and the gravitational forces, assuming that 
the gravitational interaction is only two-body. The effects of heating and ablation of the planetesimal 
as it passed through the protoplanetary envelope were included. We also allowed for the possibility of 
planetesimal fragmentation. This will occur when the pressure gradient across the planetesimal exceeds the 
material strength of the body, and the planetesimal is small enough so that self-gravity cannot counteract the 
disruptive effect of the pressure gradient (Pollack et al. 1979). Further details of the code for this calculation 
are described in Podolak et al. (1988). 

We computed trajectories for planetesimals with successively larger impact parameters with the proto- 
planet in order to determine the largest impact parameter for which planetesimal capture was possible. For 
this trajectory the planetesimal's closest distance of approach to the center of the protoplanet was defined as 
the effective radius for capture Rgf / . Of course Ref / depends on the (instantaneous) protoplanetary density 
distribution, as well as on the size, composition and random velocity of the planetesimals. 

The planetesimals were assumed to be composed of a mixture of ice, silicates and CHON with a density 
of 2.8 g cm~"^. In what follows, this mixture of ice, rock and CHON will be referred as "ice -t- rock". 
For comparison we also ran some cases with pure rock planetesimals with a density of 3.4 g cm~^. The 
planetesimal material properties are given in Table 2. 



Physical parameter 


Ice 


Rock 


M (g/mole) 
p (g cm~^) 

LogioPvapor (dynes cm~^) 


18 
1 


50 
3.4 

-24605 + 13.176 



Table 2: Properties of the planetesimal material 



The effective radius for planetesimal capture will depend not only on the size and composition of the 
planetesimals, but also on their random velocity. Following Pollack et al. (1996) we have taken the random 
velocity of the planetesimal far from the protoplanet to be 1 km s"^ in our baseline model. We have also 
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considered planetesimals with random velocities of 2 km s^^. This is aheady 15% of the keplerian velocity 
at Jupiter's heliocentric distance, and it is unlikely that the random velocities can be stirred up much higher 
than this in the short time before most of the planetesimals are captured. Fig. 3 shows the ratio of the 
effective radius for planetesimal capture to the actual radius of the protoplanet for planetesimals with radii 
of 1, 10, and 100 km composed of ice + rock and for the two values of the random velocity we considered. As 
can be seen from the figure, for the higher random velocity the effective radii decrease since it is harder to 
slow down and capture the planetesimals (especially those with larger radii). Similarly, planetesimals with 
lower velocities than in our baseline case will be captured more efficiently. 




Figure 3: Ratio of the effective radius of the planet to the actual radius as a function of time for ice + rock 
planetesimals with radii of 1 km (solid curves), 10 km (dashed curves), and 100 km (dash-dot curves). For 
each case two curves are shown. The upper curve is for a random velocity of 1 km s^^ and the lower curve 
is for a random velocity of 2 km s^"'^. 

The additional heating supplied by the influx of planetesimals is not expected to affect the course of the 
planet's evolution significantly. The radius of the protoplanet should increase as a result, and its energy 

8 



will increase as well. The temperature (both central and effective) will decrease in accordance with the well 
known negative heat capacity of a self gravitating object. The decrease of the effective temperature Te is 
more than offset by the increase of R, so that the luminosity L = AirR^asBTg increases. As a result, the 
Kelvin-Helmholz time, which is inversely proportional to the product RL will actually be somewhat shorter. 
Detailed simulations confirm this simple picture: Initially the radius for the case with planetesimal heating 
was 30% larger than for the case without planetesimal heating. After some yr the radius was only 15% 
larger, and by 10^ yr the difference was 3%. The effect of the planetesimals accreted during the evolution 
on the protoplanetary mass was negligible. 

4 Mass Accretion 

As can be seen from fig. 3, -Re// is initially much smaller than the actual physical radius of the protoplanet. 
This is because the outer layers of the protoplanet have such low density that the planetesimal cannot 
be sufficiently slowed until it enters deeper into the protoplanetary atmosphere. As a result, only those 
planetesimals hitting the protoplanet with relatively small impact parameters will be captured. 

As the protoplanet contracts, however, the density throughout most of its volume increases and the drag 
force becomes more important. Initially the effective radius increases, but once the interior densities are 
high enough, and the protoplanet continues to contract, the effective radius decreases again. The point at 
which the gas density becomes high enough to efficiently capture planetesimals depends on their size, and 
as can be seen from the figure, the turn around point is later for larger planetesimals. In addition, the 
larger planetesimals cannot be captured at all initially. Planetesimals with radii of 10 and 100 km cannot 
be captured during the first ~ 1 x 10^ and ~ 7.5 x 10^ years, respectively. Although this means that the 
time for planetesimal capture is reduced significantly for these cases, we will show that the time remaining 
is still sufficient to accrete most of the solids in the feeding zone. 

There is a complication due to the fact that encounters between planetesimals and the protoplanet 
actually take place in the presence of the Sun, and should really be treated as three-body encounters. In this 
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case, the two-body cross section needs to be multiplied by a gravitational enhancement factor: S = FgirRlj^^ 
which corrects for the three-body effects. Greenzweig and Lissauer (1990; 1992) found expressions for Fg 
based on fits to numerical simulations of the planetesimal orbits. For the conditions studied here, their fits 
give an Fg which is typically several hundred, and can be as large as 10^. However, the fits of Greenzweig 
and Lissauer were done for protoplanets small compared to the Hill sphere radius. In our case, the clump 
has an initial radius which is essentially equal to the Hill sphere radius. As a result, it is not clear that the 
numerical results of Greenzweig and Lissauer apply in this case. 

An alternative is to use the two-body gravitational cross section, nb"^ where b is the largest impact 
parameter for which a planetesimal is captured, and is computed using the planetesimal trajectory code 
described above. This however, neglects any influence due to the Sun. Fig. 4 shows the difference between 
Reff and b for different cases. As can be seen from the figure, the difference is, for the most part, a factor of 
a few. We will here use the most conservative estimate of the cross section, and take it equal to S = i^R^ff ■ 
This will overestimate the time required to accrete the planetesimals in the protoplanet's vicinity. 
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0.5 1 1.5 2 2.5 3 

Time [yrs] ^'^ 

Figure 4: Comparison of the effective radius (curves labeled "Re") and the impact parameter (curves labeled 
"b") for different planetesimal sizes as a function of time. 

The total mass of solids in the planet's feeding zone is given by 

Md^sk = Halut - a^Ht) (8) 

where Uout and a™ are the outer and inner radii of the feeding zone, respectively. Again, following Pollack 
et al. (1996) we assume cr = 10 g cm^^ in the region of Jupiter's formation, and that the planetesimals 
are spread uniformly on either side of the orbit to a distance a/, which depends on the eccentricity and 
inclination of the planetesimal orbits, as well as the Hill sphere radius of the protoplanet (see Pollack et al. 
for details). For conditions similar to those assumed by Pollack et al. a/ = 2.2 x 10^'^ cm, Uout = a + af and 
ain — a — af. The mass of solids in Jupiter's feeding zone using these values is 2.2 x 10^^ g ~ SSM^. This is 
roughly twice the estimated amount of heavy elements in Jupiter (Saumon and Guillot 2004), so that even 
if some of the planetesimals are scattered out of Jupiter's feeding zone, the protoplanet should still manage 
to capture enough material to match theoretical models of Jupiter. 
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We can now use eq. 1 to compute the rate of mass capture by the protoplanet. We consider three sizes 
of planetesimal: 1, 10, and 100 km. Since a clump is expected to start forming relatively quickly, it is likely 
that the planetesimals were in the small end of this range. Fig. 5 shows the mass of accreted material as 
a function of time for the three planetesimal sizes we considered. In all cases the planetsimals are assumed 
to be composed of a mixture of rock and ice in solar proportions. For the purpose of computing the drag 
force and the resultant ablation of the planetesimal, we view it as composed of rocky material embedded in 
a matrix of ice. When the temperature is high enough for the ice to ablate, a proportionate amount of rock 
is lost together with it. For such planetesimals, the heating and mass loss are similar to those for pure ice, 
and they are easily captured. 

jj ^ q29 Accreted Mass - Ice+Rock planetesimals 

2.5 1 \ \ \ I I \ 

1km 



1.5 ■ 



0.5 



1 0km / 

/ 

/ 

/ 

/ 

/ 

100km / 
/ 

/ 

/ 

/ 

/ 

/ 

/ 
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Figure 5: Captured mass for the three sizes of planetesimals made of ice and rock vs. time. 

1 km planetesimals can be captured almost immediately, while 10 km planetesimals can only be captured 
after the protoplanet has contracted for some tens of thousands of years. 100 km planetesimals require ~ 10^ 
yr of evolution before the average density of the protoplanet becomes high enough to capture them. In the 
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first two cases all of the available planetesimal mass is captured. In the case of 100 km planetesimals some 
90% of the available mass is captured. 

Fig. 6 shows the same calculation for the assumption of pure rock planetesimals. As expected, the 
protoplanet must contract more than for the equivalent case for planetsimals composed of a mixture of ice 
and rock, but the difference is not large. For both 1 and 10 km planetsimals all of the available mass is 
captured within ~ 10^ yr. For the 100 km planetsimals ~ 80% of the mass is captured before the protoplanet 
enters the rapid contraction phase. Since it is likely that the planetesimals will have some distribution is 
size, some of the mass will be in smaller bodies, and the fraction of solid material that will be captured 
will be even closer to 100%. Thus, assuming that none of the material is gravitationally scattered out of 
the feeding zone, a contracting protoplanet is expected to accrete most of the mass in its feeding zone well 
before it evolves to the rapid contraction stage. 

^ ^ q29 Accreted mass - Rock planetesimals 

2.5 1 \ \ \ I I \ I 



2 - / 1km 



10km 



y 



0.5 - 



0.5 1 1.5 2 2.5 3 

Time [yrs] " ' " 



Figure 6: Captured mass for the three sizes of planetesimals made of pure rock vs. time. 
The total mass of solids in the feeding zone is a direct consequence of the assumed solid surface density. 
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Using higher values of a increases the total mass of solids in the accretion zone. As it is usually presented, 
the disk instability model requires a greater surface density (by ~ 150%) for an instability to form. We 
therefore repeated our calculation using ct as 15 g cm~^. In this case, the total mass of the planetesimals in 
the feeding zone is ~ 3.3 x 10^® g 55Mearth and, as in the precvious case, essentially all the planetesimals 
are captured. We summarize our results in table 3. Here the "capture time" is the time required to capture 
essentially all of the planetesimals in the feeding zone. If this time is longer than TcoI = 3.04 x 10^ yr, the 
time for the protoplanet to reach the rapid collapse stage, then TcoI is given in the table as the capture time. 
As can be seen, all 1 km and 10 km planetesimals can be captured within a relatively short time. Things are 
slightly different with the 100 km planetesimals: here the accreted mass is ~ 95% of the solids for a random 
velocity, 14o of 10^ cm s~^. For V^o = 2 x 10^ cm s~^, less than half of the material is captured. 





cr = 10 g/cm^, Voo = 1 km/s 


0- = 10 g/cm^, Voo = 2 km/s 


cr = 15 g/cm^, V^oo = 1 km/s 


rp (km) 


Time (years) 


Mcaptured{g) 


Time (years) 


Mcapturedis) 


Time (years) 


Mcapturedig) 


1 


2.09 X 10^ 


2.19 X 102» 


4.89 X 10'' 


2.19 X 1029 


3.66 X 10^ 


3.29 X 10^9 


10 


1.52 X 10* 


2.19 X 10^9 


2.99 X 10^ 


2.19 X 10^9 


1.63 X lO'' 


3.29 X 10^9 


100 


3.04 X 10^ 


2.09 X 10^'' 


3.04 X 10^ 


9.88 X 10^^ 


3.04 X 10^ 


3.13 X 10^9 



Table 3: The capture time and accreted mass for different cases. All cases correspond to ice+rock planetes- 
imals 



5 Summary and Conclusions 

We present an evolutionary model of an isolated giant gaseous protoplanet such as might be formed from 
a clump caused by an instability in the solar nebula. Assuming that such a clump survives and becomes a 
protoplanet, we calculate the rate at which it would capture planetesimals. We considered three planetesimal 
sizes and two types of materials, ice + rock, and pure rock. The total mass of solids that is captured depends 
on the size distribution of bodies in the accretion zone and on their composition. 100 km-sizcd planetesimals 
can not be captured during the first ^ 10^ years of the evolution. This is due to the fact that until that 
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time the density and temperature throughout most of the volume of the protoplanet are too low. Bodies 
with sizes of 1 km or less can be accreted starting at times shortly after the clump's creation. Icy bodies are 
captured more easily than pure rock planetesimals. 

In the case of higher random velocities, the capture efficiency decreases and the time for capturing the 
same amount of solids is longer, but even with random velocities as high as 15% of the keplerian velocity, 
more than 16 of heavy elements can be captured. For a higher surface density (bigger amount of solids 
in Jupiter's feeding zone), the total mass of heavy elements that is captured increases proportionately. As 
can be seen from table 3 the accreted mass varies by a factor of two when the random velocity of the 100 
km planetsimals is increased from 1 to 2 km s~^. Since the latter is already 20% of Jupiter's orbital speed, 
it is unlikely that the random velocities will be significantly higher than this, especially since most of the 
planetesimals are captured before they can be highly stirred. Even in this case some 50% of the mass in the 
feeding zone is captured. If, as is likely, there is a distribution of planetesimal sizes, with 100 km being at 
the high end, then the fraction of feeding zone mass that is captured will be even higher. 

It has been suggested (e.g. Weissman 1986) that planetesimals are actually loosely bound "rubble piles". 
If this is indeed the case, they would be held together even more weakly than our ice + rock planetesimals. 
Such planetesimals would be captured more easily that those composed of ice + rock, and again, the values 
we have presented for the mass captured would be a lower limit. 

Of course all this assumes that the planetesimals remain in the feeding zone of the protoplanet and are 
not scattered out (see, e.g., Hahn and Malhotra 1999). It should be noted, however, that simulations of plan- 
etesimal scattering treat the planets as point masses, and do not include the gas drag on the planetesimal as 
it goes through the protoplanet 's extended envelope. The energy loss which results from the gas drag will act 
to inhibit scattering. In addition, the larger cross section for capture that an extended protoplanet presents, 
will also significantly increase the ratio of captured to ejected planetesimals. Thus, in the simulations of 
Hahn and Malhotra, for example, significant mass ejection requires times much longer than the 3 x 10^ years 
required for the protojupiter to accrete most of the mass in its feeding zone. Therefore, if giant planets 
are created by the disk instability mechanism, there is a good chance that many of the planetsimals will be 
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accreted before they are scattered out of the feeding zone. As a consequence, these protoplanets can collect 
a substantial amount of solids during their evolution and end up with a significantly non-solar composition. 
This finding removes an important objection to the disk instability model as a mechanism for giant planet 
formation. 

Our model assumes a non-rotating clump, in reality, the clump is likely to have some angular momentum. 
This might act to slow the contraction of the clump in its later stages, and lead to the formation of a subdisk. 
Planetesimals in this subdisk could then provide the material for forming both regular and irregular satellites. 
It is also interesting to speculate on the possibility of observational verification of our results. Some models 
of extrasolar planets whose densities are known seem to indicate non-solar compositions (e.g., Bodenheimer 
et al. 2001). The extent of the enhacement of the non-solar component may provide limits on the proposed 
accretion mechanism. However, both the fits to the observations and the models we have computed, have, 
as yet, too many free parameters, to allow useful constraints. 

There are a number of issues which may affect the conclusions presented above. In the first place, the 
captured planetesimals will contribute high-Z material to the gaseous envelope, and change its composition. 
This will also affect the opacity of the envelope. Since much of this additional material will end up in 
the convective region of the planet, this change in opacity should not have a large effect on its subsequent 
evolution. In addition, there will be some small effect due to the additional mass of the planetesimals 
themselves. Finally, the high-Z material which remains in the envelope in the form of grains will sediment 
towards the center of the protoplanet leading to the formation of a core, and this too should affect the final 
evolution of the body. We hope to address these issues in more detail in future work. 
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